#### Calculate planting elasticity

## Setup environment

setwd("~/research/coffee-dataverse/src")

do.origspec <- T
source("load.R", encoding="iso-8859-1")

output.suffix <- "-eq24rational-orig" #"-eq24-orig"

library(dplyr)
library(rstan)

## Determine production portion that is Arabica

df.variety <- df %>% group_by(Munic�pio) %>% summarize(arabica.portion=mean(ifelse(is.na(arabica.produce), ifelse(is.na(robusta.produce), NA, 0), arabica.produce) / total.produce, na.rm=T), total.average=mean(total.produce, na.rm=T))

## Calculate h-bar / b-bar

sumstats <- read.csv(paste0("sumstats", output.suffix, ".csv"))
names(sumstats)[1] <- "fullparam"
sumstats$param <- gsub("\\[\\d+\\]", "", gsub("_coeff\\[(\\d+)\\]", "_coeff\\1", sumstats$fullparam))

sumstats3 <- sumstats %>% filter(param == 'soft_harvest_part') %>% group_by(fullparam, param, region) %>%
    summarize(mean=sum(mean) / length(mean), se_mean=sum(se_mean) / length(se_mean))

df2 <- df %>% group_by(Munic�pio) %>% arrange(year, .by_group=T) %>% summarize(year=year, total.produce=total.produce, total.harvest=total.harvest, edd1000=edd1000)

## Calculate harvest portion
df2$harvest.part.bayes <- NA
for (region in unique(sumstats3$region)) {
    rows.ss <- which(sumstats3$region == region & sumstats3$param == 'soft_harvest_part')
    rows.df <- which(df2$Munic�pio == region & !is.na(df2$edd1000))
    if (length(rows.ss) != length(rows.df))
        print(paste("Skip", region))
    if (length(rows.ss) == 0)
        next

    df2$harvest.part.bayes[rows.df] = sumstats3$mean[rows.ss]
}

df3 <- df2 %>% left_join(df.variety)

prices$year.p1 <- prices$year + 1
df4 <- df3 %>% left_join(prices, by=c('year'='year.p1')) # 2000 gets 1999

df.hbp <- df4 %>% group_by(Munic�pio) %>% summarize(hpart=mean(harvest.part.bayes[!is.na(total.harvest)], na.rm=T), pavg=mean(price.arabica[!is.na(total.harvest) & !is.na(edd1000)] * arabica.portion[!is.na(total.harvest) & !is.na(edd1000)] + price.robusta[!is.na(total.harvest) & !is.na(edd1000)] * (1 - arabica.portion[!is.na(total.harvest) & !is.na(edd1000)]), na.rm=T))

load(paste0("../data/combined-betagamma", output.suffix, ".RData"))

param <- get.fit.params(output.suffix)

pp <- which(param == "cost_to_planting")

## elast = (db/bbar) / (dp/pbar) = (pbar / bbar) hbar phi = pbar hpart phi

elasts <- la$gamma[, 1, pp] * mean(df.hbp$pavg, na.rm=T) * mean(df.hbp$hpart, na.rm=T)

mean(elasts)
quantile(elasts, c(.025, .975))
